HDL lipidome quantification (from untargeted)

source("global.R")
library(ggrepel)
library(ComplexHeatmap)

lipid <- readRDS("../data/untargeted_singletp_fin4_all.rds")

Analysis based on the targeted lipidomic results. File: mx 733759_Agus_lipidomics_human plasma lipoproteins_07-2023 submit.xlsx. The peak heights were converted to ng/mL based on the amount and peak height of iSTD.

Data processing

  • Peak heights across samples were normalized using all iSTDs.

  • Features has been filtered by WCMC regarding the peak heights in samples and blanks.

  • Features with InChI Keys that couldn’t be found in the PubChem were excluded.

  • Compounds with identical InChI Keys were selected by lowest CV in pooled samples.

  • Assign lipid class from LIPID MAPS and manually curated table.

  • Concentration were calculated based the istd with most similar unsaturation and then chain length.

  • Cholesterol istd is missing. Import the value from mx733759_Agus_Lipids_Single-point quant_with cholesterol_Submit_10-17-23.xlsx

Total peak heights consists of major lipid species which can be quantified with iSTDs and minor species which don’t have corresponding iSTDs. The percentage of both types of species are shown below.

peak <- lipid$peak
lipid <- lipid$quant
peak$fdata$major <- ifelse(featureNames(peak) %in% featureNames(lipid), "Major", "Minor")
lipid.major <- summarize_feature(peak, "major")$edata
lipid.major <- apply(lipid.major, 2, function(x)x/sum(x)*100)
df <- as.data.frame(lipid.major) %>%
        rownames_to_column("Type") %>%
        pivot_longer(cols = contains("_0"), names_to = "fraction", values_to = "Perc") %>%
        mutate(fraction = factor(fraction, levels = colnames(lipid.major))) %>%
        group_by(fraction) %>%
        mutate(ypos = 100 - cumsum(Perc) + 0.5 * Perc)
ggplot(df, aes(x = "", y = Perc, fill = Type)) +
        geom_bar(stat = "identity", color = "white") +
        geom_label(aes(y = ypos, label = round(Perc, 1)), fill = "white") +
        coord_polar("y", start = 0) +
        theme_void() +
        scale_fill_npg() +
        facet_wrap(~fraction, nrow = 3)

The quantified data only included major lipid species.

CVs of pooled samples

The CVs of pooled samples for each lipid species were calculated using the following formula:

\[ CV_{lipid\ species} = SD_{lipid\ species}/Mean_{lipid\ species} * 100\% \]

We are less confident with the lipid species with a high CV. These lipid species will be excluded for further analysis.

hist(lipid$fdata$qc_cv, main = "The distribution of CVs")

sum(lipid$fdata$qc_cv>30)
[1] 11

We excluded lipid species with a CV > 30%. 11 features were excluded.

lipid <- subset_features(lipid, lipid$fdata$qc_cv<30)

The hierarchical relationship of lipid classes.

Following analysis included all quantified lipid species.

Included all quantified lipid species. 456 lipid species for analysis.

Lipid class

Pie chart

lipid.class3 <- subset_features(lipid, !is.na(lipid$fdata$
                                                       class))
lipid.class3 <- summarize_feature(lipid.class3, "class")
edata.prop3 <- apply(lipid.class3$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class3 <- lapply(sampleNames(lipid), plotPie, edata.prop3)
legend.class3 <- get_legend(pies.class3[[1]]+theme(legend.position = "bottom"))
pies.class3 <- lapply(pies.class3, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class3, nrow = 3)

Bar plot

# bar plot
plotBar(edata.prop3)

PCA

edata.scaled3 <- scale(t(edata.prop3), center = T)
res.pca <- PCA(edata.scaled3, scale.unit = T, graph = FALSE)
# fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_ind(res.pca, 
             axes = c(1,2),
             repel = TRUE,
             title = NULL
) +
        theme_cynthia_bw() +
        theme(
                legend.position = "none"
        ) +
        labs(title = NULL)

Lipid main class

Pie chart

lipid.class2 <- subset_features(lipid, !is.na(lipid$fdata$Main.class))
lipid.class2 <- summarize_feature(lipid.class2, "Main.class")
edata.prop2 <- apply(lipid.class2$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class2 <- lapply(sampleNames(lipid), plotPie, edata.prop2)
legend.class2 <- get_legend(pies.class2[[1]])
pies.class2 <- lapply(pies.class2, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class2, nrow = 3)

Bar plot

# bar plot
plotBar(edata.prop2)

Lipid category

Pie chart

# edata
lipid.class <- subset_features(lipid, !is.na(lipid$fdata$Categories))
lipid.class <- summarize_feature(lipid.class, "Categories")
edata.prop <- apply(lipid.class$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class1 <- lapply(sampleNames(lipid), plotPie, edata.prop)
legend.class1 <- get_legend(pies.class1[[1]])
pies.class1 <- lapply(pies.class1, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class1, nrow = 3) 

Bar plot

# bar plot
plotBar(edata.prop)

Following analysis only included lipids curated in LIPID MAPS.

220 lipid species for the analysis.

LM: sub class

Pie chart

lipid.class.sub <- subset_features(lipid, !is.na(lipid$fdata$LM_sub_class))
lipid.class.sub <- summarize_feature(lipid.class.sub, "LM_sub_class")
edata.prop.sub <- apply(lipid.class.sub$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class.sub <- lapply(sampleNames(lipid), plotPie, edata.prop.sub)
legend.class.sub <- get_legend(pies.class.sub[[1]]+theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 5)))
pies.class.sub <- lapply(pies.class.sub, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.sub, nrow = 3) %>%
        plot_grid(legend.class.sub, nrow = 2, rel_heights = c(3,1))

Bar plot

# bar plot
plotBar(edata.prop.sub)

PCA plot

edata.scaled.sub <- scale(t(edata.prop.sub), center = T)
res.pca <- PCA(edata.scaled.sub, scale.unit = T, graph = FALSE)
# fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_ind(res.pca, 
             axes = c(1,2),
             repel = TRUE,
             title = NULL
) +
        theme_cynthia_bw() +
        theme(
                legend.position = "none"
        ) +
        labs(title = NULL)

LM: main class

Pie chart

lipid.class.main <- subset_features(lipid, !is.na(lipid$fdata$
                                                       LM_main_class))
lipid.class.main <- summarize_feature(lipid.class.main, "LM_main_class")
edata.prop.main <- apply(lipid.class.main$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class.main <- lapply(sampleNames(lipid), plotPie, edata.prop.main)
legend.class.main <- get_legend(pies.class.main[[1]]+theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3)))
pies.class.main <- lapply(pies.class.main, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.main, nrow = 3) %>%
        plot_grid(legend.class.main, nrow = 2, rel_heights = c(4,1))

Bar plot

# bar plot
plotBar(edata.prop.main)

LM: core

Pie chart

lipid.class.core <- subset_features(lipid, !is.na(lipid$fdata$LM_core))
lipid.class.core <- summarize_feature(lipid.class.core, "LM_core")
edata.prop.core <- apply(lipid.class.core$edata, 2, function(col){
        col/sum(col)*100
})

# pie chart
pies.class.core <- lapply(sampleNames(lipid), plotPie, edata.prop.core)
legend.class.core <- get_legend(pies.class.core[[1]])
pies.class.core <- lapply(pies.class.core, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.core, nrow = 3)

Bar plot

# bar plot
plotBar(edata.prop.core)

Lipid species enriched in fractions

The analysis was based on peak heights. We assumed that the sum of peak heights in each sample represented the total lipid amount. Each lipid species was normalized to the total lipid amount before comparing across samples.

CVs of pooled samples

hist(peak$fdata$qc_cv, main = "The distribution of CVs")

sum(peak$fdata$qc_cv>30)
[1] 18

We excluded lipid species with a CV > 30%. 18 features were excluded.

peak <- subset_features(peak, peak$fdata$qc_cv<30)
major.classes <- c("Cer", "CAR", "CE", "Cholesterol", "DG", "FA", "LPC", "LPE", "PC", "PE", "SM", "TG")
classes <- peak$fdata$class
classes[!classes %in% major.classes] <- "Other"
peak$fdata$class <- classes
peak.rel <- apply(peak$edata, 2, function(x){
        x/sum(x)
})

Species level

Class level

peak.class <- summarize_feature(peak, "class")
peak.class.rel <- apply(peak.class$edata, 2, function(x){
        x/sum(x)
})
p.class <- enrich_all_fraction(bg.fraction, peak.class.rel, peak.class, TRUE)

# Heatmap of enrich species in each fraction
mat.class <- as.matrix(p.class)
mat.class[mat.class<0.01] <- 0.01
mat.class <- -log(mat.class)

cols <- RColorBrewer::brewer.pal(name = "Paired", n = 12)
cols <- colorRampPalette(cols)(13)
cols <- setNames(cols, rownames(p.class))

Heatmap(mat.class, name = "-log(p)",
        right_annotation = rowAnnotation(class = rownames(p.class), col = list(class = cols)),
        show_row_names = FALSE)